#include <hw3_common.h>
#include <cublas_v2.h>
#include <stdio.h>
#include <stdlib.h>
#include "symMul.h"

#include "symMulKernel_v1.cu"
#include "symMulCoalesceKernel.cu"

void onesMatrix(float * M)
{
    for (unsigned int i=0; i< SIZE; i++){
        M[i]=1;
    }

}

void symMulCPU(float* M, float* B)
{
    for (unsigned int i =0; i< N; i++){
        for (unsigned int j=0; j< N; j++){

            B[idx(i,j)] = 0;
            for (unsigned int k=0; k<N; k++){
                B[idx(i,j)] += M[idx(i,k)]*M[idx(j,k)];
            }

        }
    }

}

void symMul(const float* M, float * B)
{
    float * d_M;
    float * d_B;

    size_t sizeM = SIZE*sizeof(float);
    cudaMalloc((void**)&d_M, sizeM);
    cudaMemcpy(d_M, M, sizeM, cudaMemcpyHostToDevice);

    cudaMalloc((void**)&d_B, sizeM);
    //cudaMemcpy(d_B, B, sizeM, cudaMemcpyHostToDevice);

    dim3 dimBLOCK(BLOCK_SIZE, BLOCK_SIZE);
    dim3 dimGRID(GRID_SIZE, GRID_SIZE);

    symMulKernel<<<dimGRID, dimBLOCK>>>(d_M, d_B);

    cudaMemcpy(B, d_B, sizeM, cudaMemcpyDeviceToHost);

    cudaFree(d_M);
    cudaFree(d_B);
}

void symMulCoalesce(const float* M, float * B)
{
    float * d_M;
    float * d_B;

    size_t sizeM = SIZE*sizeof(float);
    cudaMalloc((void**)&d_M, sizeM);
    cudaMemcpy(d_M, M, sizeM, cudaMemcpyHostToDevice);

    cudaMalloc((void**)&d_B, sizeM);
    //cudaMemcpy(d_B, B, sizeM, cudaMemcpyHostToDevice);

    dim3 dimBLOCK(BLOCK_SIZE, BLOCK_SIZE);
    dim3 dimGRID(GRID_SIZE, GRID_SIZE);

    symMulCoalesceKernel<<<dimGRID, dimBLOCK>>>(d_M, d_B);

    cudaMemcpy(B, d_B, sizeM, cudaMemcpyDeviceToHost);

    cudaFree(d_M);
    cudaFree(d_B);
}

void symMulcublas(const float* M, float* B)
{
    float* d_M;
    float* d_B;

    size_t sizeM = SIZE*sizeof(float);
    cudaMalloc((void**)&d_M, sizeM);
    cudaMemcpy(d_M, M, sizeM, cudaMemcpyHostToDevice);

    cudaMalloc((void**)&d_B, sizeM);

    dim3 dimBLOCK(BLOCK_SIZE, BLOCK_SIZE);
    dim3 dimGRID(GRID_SIZE, GRID_SIZE);

    cublasHandle_t t_Handle;
    cublasStatus_t status;
    status = cublasCreate(&t_Handle);

    const float alpha = 1;
    const float beta = 0;

    status = cublasSsyrk(t_Handle,CUBLAS_FILL_MODE_UPPER, CUBLAS_OP_N, N, N, &alpha, d_M, N, &beta, d_B, N );

    cudaMemcpy(B, d_B, sizeM, cudaMemcpyDeviceToHost);

    cublasDestroy(t_Handle);
    cudaFree(d_M);
    cudaFree(d_B);

}


int main()
{
    float * M=(float*)malloc(SIZE*sizeof(float));
    float * Bcpu=(float*)malloc(SIZE*sizeof(float));
    float * Bgpu=(float*)malloc(SIZE*sizeof(float));
    float * Bcoalesce = (float*)malloc(SIZE*sizeof(float));
    float * Bcublas=(float*)malloc(SIZE*sizeof(float));

    genData(M, N) ;

    //make it symmetric
    for (unsigned int i=0; i<N; i++){
        for (unsigned int j=i+1; j<N; j++){
            M[idx(j,i)]=M[idx(i,j)];
        }
    }

    onesMatrix(M);

    genData(Bcpu, N);
    genData(Bgpu, N);
    genData(Bcublas, N);
    genData(Bcoalesce, N);

    double begin = wtime();
    // matrix multiply with cpu
    for (unsigned int i=0; i< ITER; i++){
        symMulCPU(M, Bcpu);
    }
    double end = wtime();
    double cputime = (end-begin)/(double)ITER;

    // matrix multiply with kernel
    begin = wtime();
    for (unsigned int i=0; i< ITER; i++){
        symMul(M, Bgpu);
    }
    end = wtime();
    double gputime = (end-begin)/(double)ITER;

    // matrix multiply with coalesce
    begin = wtime();
    for (unsigned int i=0; i< ITER; i++){
        symMulCoalesce(M, Bcoalesce);
    }
    end = wtime();
    double coalescetime = (end-begin)/(double)ITER;

    //matrix multiply with cublas
    begin = wtime();
    for (unsigned int i=0; i<ITER; i++){
        symMulcublas(M, Bcublas);
    }
    end = wtime();
    double cublastime = (end-begin)/(double)ITER;

    float max1 = 0;
    float max2 = 0;
    float max3 = 0;
    float max4 = 0 ;
    for (unsigned int i=0; i < N; i++){
        for (unsigned int j=0; j< i+1; j++){
            float diff1 = abs(Bcpu[idx(i,j)]-Bgpu[idx(i,j)]);
            float diff2 = abs(Bcpu[idx(i,j)]-Bcublas[idx(i,j)]);
            float diff3 = abs(Bgpu[idx(i,j)]-Bcublas[idx(i,j)]);
	    float diff4 = abs(Bcpu[idx(i, j)] - Bcoalesce[idx(i, j)]);
            if (diff1 > max1){
                max1 = diff1;
            }
            if (diff2 > max2){
                max2 = diff2;
            }
            if (diff3 > max3){
                max3 = diff3;
            }
	    if (diff4 > max4){
		max4 = diff4 ;
	    }
//           printf("row:%d col:%d cpuvalue:%f gpuvalue:%f cublasvalue:%f\n", i, j, Bcpu[idx(i,j)],Bgpu[idx(i,j)], Bcublas[idx(i,j)], );
        }
    }

    printf("cpu-gpu:%f  cpu-cublas:%f  gpu-cublas:%f    coalesce:%f\n", max1, max2, max3, max4);
    printf("cpu time:%f  gpu time:%f   cublas time:%f    coalescetime:%f\n", cputime, gputime, cublastime, coalescetime);


    free(M);
    free(Bcpu);
    free(Bgpu);
    free(Bcublas);

}
